
Criteria for Best Practice

Witch Survey
Analysis

library(terra)
# Calculate distance to forest
forest <- terra::segregate(nlcd, classes = 42) # Extract the forest layer
forest <- terra::classify(forest,
rcl = matrix(c(1, 0, 1, NA),
nrow = 2,
ncol = 2)) # Reclassify 0 as NA
forest <- terra::distance(forest)
forest <- project(forest, "EPSG: 4326") # Reproject to WGS84
forest <- mask(crop(forest, ext(tetlin)), tetlin)
names(forest) <- "forest"
# Calculate distance to water
water <- terra::segregate(nlcd, classes = 11) # Extract the water layer
water <- terra::classify(water, rcl = matrix(c(1, 0, 1, NA),
nrow = 2,
ncol = 2)) # Reclassify 0 as NA
water <- terra::distance(water)
water <- project(water, "EPSG: 4326") # Reproject to WGS84
water <- mask(crop(water, ext(tetlin)), tetlin)
names(water) <- "water"

# Required packages
library(sf)
library(terra)
# Get witch observation data from ServCat
sites <- pull_witch_data() %>% filter(year == 2024)
# Reformat for `unmarked` package
sites <- reformat(sites)
# Extract covariates to sites
sites <- data.frame(sites,
forest = terra::extract(forest, sites)$forest,
water = terra::extract(water, sites)$water)# Required package
library(unmarked)
# Create an unmarked data frame (scaled covariates)
unmarked_df <- unmarkedFrameOccu(y = sites[,5:13], # observations
siteCovs = sites[,4:5]) # covariates
sc <- scale(site_covs) # scale them
siteCovs(unmarked_df) <- sc # add back
# Fit single-season occupancy model
mod <- unmarked::occu(formula = ~forest ~ water +
forest,
data = unmarked_df[[1]])
# Look at estimates
mod@estimatesOccupancy:
Estimate SE z P(>|z|)
(Intercept) 0.0923 0.213 0.433 0.6651
water 0.8643 0.337 2.568 0.0102
forest -0.4465 0.300 -1.487 0.1370
Detection:
Estimate SE z P(>|z|)
(Intercept) -0.0507 0.115 -0.44 0.6596
forest 0.7406 0.325 2.28 0.0226

Criteria for Best Practice
Witch Survey
Results
# Calculate predicted occupancy
pred <- rbind(unmarked::plotEffectsData(mod, "state", "forest"),
unmarked::plotEffectsData(mod, "state", "water")) |>
mutate(covariateValue = case_when(
covariate == "forest" ~ covariateValue * attr(sc, 'scaled:scale')[[1]] + attr(sc, 'scaled:center')[[1]],
covariate == "water" ~ covariateValue * attr(sc, 'scaled:scale')[[2]] + attr(sc, 'scaled:center')[[2]]
))#' Create a leaflet map of occupancy within a refuge
#'
#' @param ras a `terra::ras` raster of psi estimates
#' @param s an `sf::st_point` of the sites surveyed
#' @param r an `sf` multipolygon of the refuge boundary
#' @para p whether to map the predicted value of psi ("Predicted") or SEs ("SE")
#' @param h the height of the leaflet map returned
#' @param w the width of the leaflet map returned
#'
#' @return a leaflet map
#'
#' @import RColorBrewer
#' @import leaflet
#' @import terra
#'
#' @export
#'
#' @example
#' \dontrun{
#' create_map(ras = psi, s = sites, r = tetlin, p = "Predicted", h = 650, w = 300)
#' }
create_map <- function(ras,
s,
r,
p = "Predicted",
h = NULL,
w = NULL) {
if (p == "Predicted") {
x <- c(round(minmax(ras)[[1,1]],2),
round(minmax(ras)[[2,1]],2))
grp <- "Psi"
ras <- ras$Predicted
} else if (p == "SE") {
x <- c(round(minmax(ras)[[1,2]],2),
round(minmax(ras)[[2,2]],2))
grp <-"SE"
ras <- ras$SE
}
pal_rev <- colorNumeric(RColorBrewer::brewer.pal(5, "Spectral"),
x,
reverse = TRUE,
na.color = "#00000000")
pal <- colorNumeric(RColorBrewer::brewer.pal(5, "Spectral"),
x)
leaflet(height = h, width = w) |>
addTiles() |>
addRasterImage(ras,
colors = pal_rev,
maxBytes = 10168580,
opacity = 0.75,
group = grp) |>
addCircleMarkers(data = s, lat = ~Y, lng = ~X,
radius = 0.5,
color = "black",
group = "sites") |>
addPolygons(data = r,
fill = FALSE,
color = "black",
group = "Tetlin",
weight = 0.5) |>
addLayersControl(overlayGroups = c(grp,
"sites",
"Tetlin"),
options = layersControlOptions(collapsed = FALSE)) |>
addLegend(pal = pal,
values = x,
title = grp,
labFormat = labelFormat(transform = function(x) sort(x, decreasing = TRUE))) |>
addMiniMap(height = 100,
width = 100) |>
addScaleBar()
}
base_map <- function(s, r, h = NULL, w = NULL) {
leaflet(height = h, width = w) |>
addTiles() |>
addCircleMarkers(data = s, lat = ~Y, lng = ~X,
radius = 0.5,
color = "black",
group = "sites") |>
addPolygons(data = r,
fill = FALSE,
color = "black",
group = "Tetlin",
weight = 0.5) |>
addLayersControl(overlayGroups = c("sites",
"Tetlin"),
options = layersControlOptions(collapsed = FALSE)) |>
addMiniMap(height = 100,
width = 100) |>
addScaleBar()
}# Source the `create_map()` function
source("./R/create_map.R")
# Import data
tetlin <- sf::st_read("data/shapefile/tetlin.shp",
quiet = TRUE)
sites <- read.csv("data/csv/sites.csv")
psi <- terra::rast("data/raster/psi/psi.tif")
# Create a map
create_map(ras = psi,
s = sites,
r = tetlin,
p = "SE",
h = 650,
w = 300)

Criteria for Best Practice



Witch Survey
Report
Quarto Code
```{r}
---
title: "Tetlin Witch Report"
author: Jane Biologist
format: html
fig-align: center
editor: source
---
```{{r setup}}
#| echo: false
#| message: false
knitr::opts_chunk$set(warning = FALSE,
echo = FALSE,
message = FALSE,
fig.retina = 3,
fig.align = "center")
library(unmarked)
library(terra)
library(tidyverse)
library(RColorBrewer)
library(sf)
library(leaflet)
```
```{{r load_data}}
#| cache: true
# Load site data
dat <- read.csv("data/csv/dat.csv")
source("R/simulate_data.R")
source("R/create_map.R")
# Scaled covariates
sc <- dat |>
dplyr::select(forest, water) |>
scale()
# Load site data and scale them
load("./data/rdata/unmarked_df.Rdata")
```
```{{r fit_model}}
#|cache: true
# Fit single season occupancy model
mod <- fit_model(unmarked_df)
```
## Introduction
Invasive witches have become a management concern at Tetlin National Wildlife Refuge. As such, there is a need to estimate witch occurrence within the Refuge.
## Methods
### Data Collection
We visited a sample randomly distributed sites across Tetlin Refuge. At each site, we spent one hour looking and listening for witches. We revisited each site eight times.
### Model
We estimated witch occupancy and detection using a single-season occupancy model. We used the `unmarked` R package. [blah, blah, blah]
## Results
We surveyed a total of `r nrow(dat)` sites. The average distance to water at our sites was `r round(mean(dat$water), 2)` m. The average distance to forest at our sites was `r round(mean(dat$forest), 2)` m.
```{{r}}
#| out-width: "50%"
#| fig-cap: "A map of the sites surveyed for witches, Tetlin National Wildlife Refuge, Alaska."
# Import data
tetlin <- sf::st_read("data/shapefile/tetlin.shp",
quiet = TRUE)
sites <- read.csv("data/csv/sites.csv")
# Create leaflet map
base_map(sites, tetlin)
```
We observed witches on `r sum(dat[6:13])` of 800 site visits, for a naive occupancy of `r round(sum(dat[6:13])/ncell(dat[6:13]), 2)`.
```{{r plot_psi}}
#| fig-height: 3
#| fig-cap: Occupancy of witches at Tetlin National Wildlife Refuge, Alaska, 2024.
# Calculate predicted values (unscaled)
pred <- rbind(plotEffectsData(mod, "state", "forest"), plotEffectsData(mod, "state", "water")) %>%
mutate(covariateValue = case_when(
covariate == "forest" ~ covariateValue * attr(sc, 'scaled:scale')[[1]] + attr(sc, 'scaled:center')[[1]],
covariate == "water" ~ covariateValue * attr(sc, 'scaled:scale')[[2]] + attr(sc, 'scaled:center')[[2]]
))
# Plot predicted values (psi)
ggplot(pred, aes(x = covariateValue, y = Predicted),
group = covariate) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
linetype = 2,
alpha = 0.1) +
xlab("Distance (m)") +
ylab("Psi") +
facet_grid(~covariate, scales = "free")
```
```Rendered Report

